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Summary 

This is a user’s manual of the CMOTT turbulence module, version 2.0, developed for 
the NPARC code. The module is written in a self-contained manner so that the user can 
use any turbulence model in the module without concern as to how it is implemented and 
solved. Three two-equation turbulence models have been built into the module: Chien, 
Shih-Lumley and CMOTT models, and all of them have both the low Reynolds number 
and wall function options. Unlike Chien’s model, both the Shih-Lumley and CMOTT 
models do not involve the dimensionless wall distance in the low Reynolds number 
approach, an advantage for separated flow calculations. The Van Driest transformation 
is used so that the wall functions can be applied to both incompressible and compressible 
flows. The manual gives the details of the turbulence models used and their numerical 
implementation. It also gives two application examples, one for subsonic and the other 
for transonic flow, for demonstration. The module can be easily linked to the NPARC 
code for practical applications. 
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1 Introduction 

This manual describes the version 2.0 of the CMOTT turbulence module developed at 
the Center for Modeling Of Turbulence and Transition (CMOTT) for the NPARC code 
(Cooper and Sirbaugh, 1989 and 1990). The present version differs from the previous 
one (Zhu and Shih, 1995) in the following two major aspects: 1) use of the delta form of 
turbulence transport equations, and 2) inclusion of the wall function approach. 

The purpose of developing the turbulence module is to e nhan ce the turbulence modeling 
capability of NPARC. The module is written in such a way that it, on the one hand, can be 
easily linked to NPARC for practical calculations, and on the other hand, can be updated 
in time to include the state of the art of turbulence models suitable for applications in 
aerospace and aero-propulsion systems. With the aid of the module, turbulence model 
developers can also take the advantage of the well-established and sophisticated CFD 
code to test turbulence models under development for a variety of complex flows which 
are intractable with simple research codes. 

Under the widely used Boussinesq’s isotropic eddy- viscosity concept, the Reynolds- 
averaged equations governing turbulent flows are of the same form as those governing 
laminar flows, except that the laminar viscosity fi is replaced by the effective viscosity 

faff — f 1 + fa (1) 

Therefore, a mean flow solver can be used to calculate turbulent flows once the turbulent 
eddy-viscosity fit is available. For most CFD codes, especially those for compressible 
flow calculations, the laminar viscosity fi is a variable not a constant. In this case, few 
changes are required for a mean flow solver to use the turbulence module. The input 
to the module are the mean flow variables, boundary and geometric information which 
are to be provided by a mean flow solver. The output of the module are the turbulent 
eddy- viscosity fi t and relevant turbulent source terms which are needed for the mean flow 
calculation. The interaction between the mean flow solver and the turbulence module will 
give the final turbulent flow solution. 

In the module, the three low Reynolds number K-e turbulence models have been im- 
plemented: Chien (1982), Shih-Lumley (1993) and CMOTT realizable models (Yang et 
al., 1995; Zhu and Shih, 1995). Chien’s model is one of the well-known low Reynolds 
number K-e models. However, it has some undesirable deficiencies. First, a near wall 
pseudo-dissipation rate is introduced to remove the singularity in the dissipation rate 
equation at the wall. The definition of the near wall pseudo-dissipation rate is quite ar- 
bitrary. Second, the model constants are different from those of the standard K-e model 
(Launder and Spalding, 1974), making the near wall model less capable of ha n dli n g flows 
containing both high Reynolds number turbulence and near wall turbulence. Patel et 
al. (1985) required as the first criterion, the ability of the near wall models to be able 
to predict turbulent free shear flows. Third, the dimensionless wall distance y + is used 
in the damping function for the eddy- viscosity. Because the y + involves the friction 
velocity U T which is equal to zero at separation or reattachment points, any model using 
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y + may have difficulties for separated flows. The Shih-Lumley and CMOTT models are 
free of these deficiencies. The two models differ from one another in the formulation, 
one using the standard constant and the other a new formulation. The new has the 
following desirable features: (a) It is derived from a rigorous realizability analysis (Shih et 
al., 1994) that requires the non-negativity of the turbulent normal stresses and Schwarz’ 
inequality between any fluctuating quantities. As a result, unlike most of the existing 
models, it satisfies the realizability conditions, (b) It accounts for the effect of the mean 
deformation rate by which the eddy- viscosity will be significantly reduced to an adequate 
level to mimic complex flow structures, (c) It is easier to use, as compared with other 
formulations. Simplicity is of great value for practical engineering applications. Successful 
applications of this new C M formulation can be found in Shih et al. (1994 and 1995). The 
module with these turbulence models has been applied to a number of flows including 
flows over a flat plate, in an ejector nozzle, in a transonic diffuser, and a boat-tail nozzle 
flow (Yang et al., 1995). For all the flow cases tested so far, it has been found that both 
the Shih-Lumley and CMOTT models produce improved or similar predictions compared 
with the Chien model. The CMOTT model with variable turns out to be more com- 
putationally robust than the other two. It was able to give numerical solutions in cases 
where the models with constant suffered from numerical instability. 

The major problems or difficulties associated with the low Reynolds number turbulence 
models are: 1) They require very fine grid spacing in near-wall regions, thus increasing 
considerably computational burden, especially in three-dimensional cases. Moreover, the 
highly stretched nature of mesh distribution may have an adverse impact on numerical 
stability. 2) Most of models are not of tensorial invariant form, that is, they contain a dis- 
tance parameter normal to the wall. The wall-distance dependency causes inconvenience 
for model applications in complicated geometries. Currently, great effort is being given in 
the area of near-wall turbulence modeling to remove this dependency, but no satisfactory 
result has been obtained yet. 3) Most of the low Reynolds number turbulence models 
were fine-tuned against attached flows, which is, of course, not sufficient to guarantee 
their good performance for separated flows. 

An alternative is to use the high Reynolds number turbulence models. Here, the 
governing equations are integrated to a point far outside the viscous sublayer rather 
than down to the wall, and the near-wall region is bridged over with the wall functions. 
Although in theory the wall function approach is only valid for certain attached flows with 
no pressure gradient and mass transfer, it has been applied in practice to many separated 
flows with varing degree of success. For those flows where maximum shear stresses occur 
far away from the wall, the near-wall turbulence modeling is not critical for overall flow 
simulations. In these cases, use of the wall functions has a very beneficial effect on the 
stability and economy of computations. Although the principle argument for originally 
adopting the wall function approach (economy of grid points) has been weakened as larger 
and faster computers have become available, it will still find its applications in predicting 
complex flows, especially for large scale engineering problems. 

In the present work, we extend the turbulence module by including the wall functions. 
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For incompressible flows, the universal law of the wall may be expressed as 



where k = 0.41 and C = 5.2. The derivation of Eq. (2) is based on the assumption that 
the shear stress in the region close to the wall is constant and equal to the wall shear stress. 
It has been shown (Viegas et al., 1985; Huang and Coakley, 1993) that the same form 
also exists for compressible flows with the velocity U being replaced by the Van Driest 
transformed velocity U c (Van Driest, 1951). For the K-e turbulence models, the convection 
and diffusion terms of their transport equations are negligible in the inertial sublayer so 
that local equilibrium prevails, which implies that the production of the turbulent kinetic 
energy K is equal to the dissipation rate e of K. The local equilibrium condition leads 
to two simple relations which can be used as boundary conditions for K and e for both 
incompressible and compressible flows. 

Although the wall functions look simple, their numerical implementation is not trivial. 
The main difficulty comes from the logarithmic law in which both U and U T are unknown, 
and U r cannot explicitly be solved for. It is prone to being numerically unstable if one uses 
Eq. (2) and iteratively solves U T to obtain the boundary conditions for the Navier-Stokes 
equations. In the module, we use an implicit procedure which directly incorporates Eq. (2) 
into the Navier Stokes equations. In this way, there is no need to solve Eq. (2) for U T 
by sub-iteration. The implicit method turns out to be more stable than the explicit one. 
Another important issue is the artificial viscosity. Chitsomboon (1994) found that the 
artificial viscosity originally implemented in the NPARC code totally spoiled the solution 
of the wall functions. This was because the artificial viscosity became unrealistically large 
in the vicinity of walls due to very steep velocity gradients resulting Horn the coarseness 
of grid spacing as required by the wall function approach. He fixed up this problem 
by extrapolating velocities at the wall rather than using the physical values of no-slip 
velocities, when calculating the artificial viscosity. In the module, we simply turn off the 
artificial viscosity in the near-wall region. 

In general, turbulence model equations require special treatment to ensure numerical 
realizability such as the positiveness of K and e. They are also often of source- dominate 
nature, which sometimes makes the linearization of source terms crucial for computational 
stability. Due to these considerations, we used the non-delta form of the transport equa- 
tions in the previous version of the module (Zhu and Shih, 1995). The non-delta form 
leads to simpler linearization and is more effective to ensure the positiveness of the turbu- 
lent kin etic energy and its dissipation rate than the delta form. However, the non-delta 
formulation requires 8 additional arrays in the three-dimensional code, and the above ad- 
vantages did not manifest themselves to be obvious in all the test cases computed so far, 
which may be due to the fact that, on the one hand, the stability of solution process is 
constrained largely by the mean flow solver, and on the other hand, positivity for models 
of K-e type can be ensured by simple numerical measure such as clipping. Therefore, the 
present module uses the delta form, which is in line with the NPARC code. 
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2 Turbulence Models 
2.1 General Form 

In accordance with the NPARC code, a non-dimensional form of equations is adopted. The 
three low Reynolds number two-equation turbulence models built into the module have 
the following common form in which the Reynolds stresses t»j(= —pu{Uj) are calculated 
by 

Til = + v iti - \v k _>) - \pK6ij ( 3 ) 

where R e is the reference Reynolds number. The turbulent eddy viscosity fi t , the turbulent 
kinetic energy K and its dissipation rate e are calculated by the following equations 

Pi = R.UC„pK 2 /e ( 4 ) 

(p K ),, + \pVjK - R;\p + = P-pt + D (5) 

(pt) f + [pU jC - RT,'(p + = hCi^P - f,c 2 pj 4 - E ( 6 ) 

where 

P = TtiUtj ( 7 ) 

h = 1. h = 1 - 0.22exp[— (JZt/6) 1 ], R, = (8) 

fie 

P is the turbulent production. For high Reynolds number models, f tl = f 1 =f 2 = l and 
D = E = 0. The differences in the three models are given below. 


2.2 Low Reynolds Number Form 

Chien’s K-e model. 


C„ = 0.09, Ci = 1.35, C 2 = 1.8, <t k = 1, <r e = 1.3 (9) 

fn = 1 exp(— 0.0115y + ), y + = ^ ef>UT ^ n (10) 

D=-2flrV* (u) 

yl 

E = _^^! eX p(-0.5 y + ) (12) 

j/n 


6 



At the wall, the values of K and e are both set to zero. 


Shih-Lumley’s K-e model. 


ai 


1.44, C 2 — 1.92, CTjf — 1, <x € — 1.3 

(13) 

Rk = heXp* 

r 

(14) 

I" 9 , as = 5 * IO - 10 


0 

II 

(15) 


(16) 

, S« = 

(17) 


The wall boundary conditions for K and e are 

K = 0.250u*, e = 0.251i?e 


PK 


(18) 


CMOTT K-e model. This model is the same as the Shih-Lumley model, except that 
Cfj. is calculated by 


where 


= min[0.09, (Ao + A,U*K/e ) ] 


A 0 = 4, A, = V 6 cos <j) 


(19) 

( 20 ) 


U' = yjS'^ + fitf ft* 

i ShShSL 

4>= -arc cos (v^W), W= 

S* = yj S*j S*j , fit j = 2 — ^i.*) 

In the above formulas, j/ n refers to the normal distance from the wall. 


( 21 ) 

( 22 ) 

(23) 
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2.3 Wall Functions 


The compressible law of the wall (Huang and Coakley, 1993) is used in the turbulence 
module. Following the NPARC nondimensionalization, this law can be written as 


* - r - 

(24) 

where 




Ur — \J{j 1 P^wall 
= RtU r y{pl n)v,au 

(25) 


k = 0.41, E = 8.4317 


and U c is the Van Driest transformed velocity defined by (Van Driest, 1951): 


'S 

II 

to* 


(26) 

where 

3 

0*1 IK 
II 



_ 2 T wa u 

~ (7 - 1 )P'< 

(27) 


D = VA 2 + B 


In the near-wall region, with the convection neglected the energy equation 
reduced to give an expression for the total heat flux 

can be 


9 = qwaii + Ur 

(28) 

and with the local equilibrium assumption (Launder and Spalding, 1974), the turbulent 
kinetic energy K and its dissipation rate e can be calculated by 


r ti7 all/ P 

0.3 

(29) 


c _ ( Twaul'pf 1 2 
ny 

(30) 


In the above expressions, the subscript wall refers to the value of the corresponding 
function at the wall. Equations (24) and (28) - (30) form the wall functions which are 
used to bridge over the first grid point and the solid wall. 
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3 Calculation Procedure 

The following presentation is only restricted to the numerical aspects related to the tur- 
bulence module. Refer to Cooper and Sirbaugh (1989 and 1990) about the details of the 
NPARC code. 

3.1 General Form of Turbulence Equations 

The turbulent transport equations ( 5) and ( 6) in general curvilinear coordinates 0 
may be written in the following form: 


dt Q + d{Fi + d v F2 + d^Fz = d^G\ + drjGz + d^Gz + S 

(31) 

Q = J- 1 

P K 

p£ 

1 F t — 1 I" pV,K 1 „ _ j f p K V(i-VK ' 

> Fi ~ 3 pU.t j ' ‘ [ p.V(i ■ Ve 

(32) 



II 

£ 

II 

£ 

fcT 

11 

£ 

(33) 



6 -L 6 = y, 6 = C 

(34) 



»■.-**<*+£) 

(35) 



p, = r; 1 {p+^) 

(36) 

V£ • V<f> 

= 

tii + « + w 

(37) 


+ 

+ £y*7y + izVz)9r,<t> + (C*Cx + CyCy + (zCz)dc4> 

V 77 • V<f> 

= 

till + V 2 y + Vl)dv<f> 

(38) 


+ 

tiixix + Vv(v + VzL)d{4> + til . *Cx + VvCy + Vz(z)dc(f> 

vc ■ v<£ 

— 

til + C y 2 + CW 

(39) 


+ 

(CxCx + CyCy + + ti*Vx + CyVy + C zVz)dr,<f> 


J is Jacobian of coordinate transformation, ?7, V and W are the contravarient velocities, 
<f> = K or e, and 5 contains the source terms in Eqs.( 5) and ( 6), respectively. 
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3.2 Time Differencing 

By using Euler backward differencing, Eq.( 31) can be written as 


A Q + A t[d(F? +1 + dr,F? +1 + d ( F? +1 - ( d + d„G^ +1 + d ( G% +1 + S n+1 )] = 0 (40) 
where 


A Q = Q n+1 - Q n 


(41) 


At = t n+1 - t n 

and the superscript n refers to the time level. 


(42) 


3.3 Time Linearization 

Eq.( 40) can be linearized by expressing all flux vectors in terms of the conservation 
variable vector Q 


pn+l = F n + AtdtF . = pn + ( 43 ) 

G? +1 = G? + A tdtGi = GT + Bi AQ (44) 


where 


S n+ 1 = S n + A tdtS = S n + CAQ 



Ui 0 ' 
0 Ui 


(45) 


(46) 


Bi = 


dGi 

dQ 


i o 

0 &22 


(47) 


dS Cll C 12 

&Q c 2 i c 22 


(48) 
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6 J 1 = [a x 3 € + J 1 hkUxVx 4 tv Vv + tzVz)dr, 4 J Vtf(£xC* 4 £ v Cy + tztz)dc\(J I p) 
b\ 2 = 4 J ^pe{txVz 4 tvVv 4 tzVz)9^ 4" J l p-c(txCx 4 tvCy 4" tztz)9(]{J/ p') 

b h = [ a2d v 4- j~* PK(Vxtx + Vvtv 4 Vztz)d( 4 J^PKiVxCx + VvCv + Vztz)d(](J/ p) 
b\ 2 = \P 2 d v + J-'fi'iVxtx 4- Vvtv + V ztz)di 4 J'V^xC* + Vvtv + Vztz)d c }(J/ p) 

&11 = [a 3 ^( 4 J~ 2 pK{txtx 4- tvtv 4 (ztz)&t + J~ X Pk{txV X 4 tvVv + C**?*)^tj](«7/ p) 
b 2 2 = [y3 3 3 C + J* 1 Pcitxtx 4 tvtv 4 tztz)d{ 4 j 2 Pc(txVx 4 tvVv 4 tzVz)Sr)\(J/ p') 


a* = +e,+ n), p = j-wn +<?+«) 

<x' = J-^Avl + ’il+’il), P-J-'r.i.vl + vl + vl) ( 50 ) 

a 3 = J-V«(G + G + Cl), = J-V.(G + e + (?) 

After introducing Eqs.( 43)-( 45), Eq.( 40) becomes 

{7 + A t[d ( A 1 + d v A 2 4 - d ( A 3 - {d t Br 4 d „£ 2 + d ( B s + C)]}AQ = RHS (51) 

where 

RHS = At(-3*2? - d v F. 2 " - d c F? 4 + d v G% + d c G” 4- S n ) (52) 

7 is an identity matrix. 


3.4 Approximate Factorization 

By using the approximate factorization method, the three-dimensional equation ( 51) can 
be reduced to the following three unidimensional equations 


[7 4- A«(d { Ai - d € 7?i)]A£?** = RHS 


[I + At(d„A 2 - dr,B 2 )\AQ * = A Q" (54) 


[7 + At(d c A 3 - d c B 3 - C))AQ = A Q* 

Solving sequentially these equations finally gives the solution of Eq.( 31) 


in+l _ n n 


= Q n + AQ 


(55) 

(56) 
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3.5 Space Differencing (LHS) 

For steady-state computations, differencing the left-hand side (LHS) of Eqs.( 53)-( 55) is 
not crucial to the accuracy of final solutions. Therefore, the first-order accurate upwinding 
is used for the inviscid flux terms and the central-differencing is used for all the other 
terms. 

With the first-order upwinding, the inviscid terms in Eq.( 53) can be written as 



A i — A + -f A~ A + — 0 a~ — U 0 

Q JJ+ 3 ^1 Q JJ- 

(57) 


P + = 0.5(17 + |t7|), U~ = O.S(<7 - |P|) 

(58) 


d t VA q' = Uf Aq} - P/^A,)., + £7- +1 A,j +1 - U-Aq) 
dfUAq’ = P/A q] - P/^A,?., + - U ;H 

(59) 

where Aq 1 and A q 2 are the components of the intermediate vector A Q**. 

By using the central- differencing and neglecting cross-derivative terms in 
viscous terms in Eq.( 53) can be written 

Eq.( 49), the 


a f 6J,A 9 ' = 0.5{(a? +1 + a)){(J / p), tl Aq) +1 - (J/p)jAqj] 



-(a] + oJ.jJK J/p)jAq} - (J/p)j-iAqj_ 1 ]} 

(60) 


d t b\,Aq' = 0.5 ■{(#„ + 0})[V/P)i» A 9 ? +1 - (J/p)jAql] 


~(/3J + 


After in 

troducing Eqs.( 59) and ( 60) into Eq.( 53), we have 



AfA,!^ + Bj'Aqj + Cj'A ,j +I = r 1 
AfAq)_ t + BpAq] + Cf Aql +1 = r J 

(61) 

where 

Af = -At[p;_, + 0.5(0?., + o})^.,] 

Bf = 1 + Ai[P/ - Pf + 0.5(a? +1 + 2a? + a?_,)(J/ri/ 
C? 1 = Ai[P/ +1 - 0.5(a? +1 + «})(Jlph» ] 

(62) 


12 



( 63 ) 


Af = -A t\Vf_ x + M{PU + #)(• 

Bf = 1 + A t[Uf - U~ + 0.5(# +1 4- 2/3j + PU){J/p)j] 

Cf = A t[Uf +1 - 0.5(# +1 + #)(J/p); + i] 

and r 1 and r 2 are the components of the vector RHS. 

The same procedure can be used for Eqs.( 54) and ( 55), except that for Eq.( 55), the 
coefficients B in Eqs. ( 62) and ( 63) are replaced by 


B l 1 = Bf 1 - c n At 
Bf 2 = Bf 2 - c u At 
Bf 1 = Bf 1 - c 21 At 
Bf 2 = Bf 2 - c 22 At 


(64) 


3.6 Space Differencing (RHS) 

To ensure both accuracy and stability of numerical solutions, the hybrid linear/parabolic 
approximation (HLPA) scheme (Zhu, 1991) is used to calculate the convection terms of the 
right-hand side (RHS) in Eq.( 52). It has been shown (Zhu, 1992) that the HLPA scheme of 
second-order accuracy works nearly as well as the third-order accurate SMART (Gaskell 
and Lau,1988) and SHARP (Leonard, 1988) schemes in terms of eliminating numerical 
diffusion while retaining boundedness of numerical solutions. Consider the inviscid flux 
in the £ direction. The HLPA scheme evaluates the value of variable <f> at the cell-face 
j — 1/2 as follows 


4>j- 1/2 — ^-1/2^'— 1 + Uj-l/ifa ^ 


(65) 


where 


8< f> = + U i- - 3 £ +1 (66) 


1 if + <t>j-2)/(<f>i ~ 4>j- 2)1 < 1 

^ -1 ^ 2 1 0 otherwise 


1 /. — 


f 1 -2<f>j + <l)j+i)/(<f>j-i - < 1 

i~ x t 2 I 0 otherwise 


(67) 


13 



and. U + and U~ are defined in Eq.( 58). 

It can be seen that Eq.( 65) is in fact the result of the first-order accurate upwinding 
with an additional term 84> added. The additional term may be viewed as an antidiffusive 
correction to the upwinding. The conventional central differencing scheme is used to 
calculate all the other terms in RES. 


3.7 Source Terms 

For the source term in Eq.( 31), we have 


S = 


Si 


51 — J~ l \P — pe + D\ 

5 2 = J~'[(C 1 P-f 2 C 2 pe)e/K + E] 


( 68 ) 


and the elements in the source term Jacobian matrix C in Eq.( 48) axe calculated by 


cn = 


dSi 

dQi' 


Cl 2 = 


dS i 
dQ 2 ’ 


C 2 1 = 


dS 2 
dQ i’ 


c 22 = 


dSi 

dQ 2 


where 


Qi = J-'pK, Q 2 = J-'pe 

From the stability requirement for the type of Eq.( 61), we have 


(69) 


(70) 


Cll , Ci 2 , C 21 , C 22 <0 (71) 

Noting that certain degree of freedom exists in dealing with the left-hand side terms of 
Eq.( 51) for steady-state calculations, we tested the following two linearization methods: 


Method 1. By ’’exactly” calculating the derivatives in Eq.( 69), we obtain 


P d(J-'D) 

11 2 pK + 8Q 1 

P d(J-'D) 
Cl2 = ___ i + __ 


4fi| 


d(J~ 1 E) 


C!1 = ClP ?x5 +C2 ^ 1/2 + "36 (/! “ 1)1+ 8Q, 

r « tf ,n lt n , , 

Cu - + 36 1 V* ~ m + 8Q, 

which do not satisfy the stability requirement Eq.( 71). 


(72) 
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Method 2. Neglecting the D and E terms, we can write the source term as 

5a = J~'P - j^Qi 

S 2 = ^ Q 2 

and the derivatives can be formally written as 


( 73 ) 


Cll = 


e 


K 


c i2 = 0 
C 2 1 = 0 


(74) 


I „ _ fiCi e 

{ c 22~ K 

which meet the stability requirement Eq.( 71). 

Numerical tests conducted so far have shown that both methods produce little differ- 
ence in terms of both convergence rate and accuracy. Therefore, only the method 2 which 
leads to a rather simple implementation is adopted in the present module. 


3.8 Wall Function Implementation 

In the NPARC code, the following nondimensional, Reynolds-averaged Navier-Stokes 
equations are solved: 

3Q , dF 1 , dF 2 dFs _ dG 1 dG 2 dGz . . 

dt + di + dT) + 8( d( + drj + d( ( ' 

where Q, Fj and Gj (j = 1,2,3) are the conservation variable vector, in viscid flux vectors 
and viscous flux vectors, respectively. Because only the viscous flux vectors need to be 
modified with the use of the wall functions, their detailed forms are given below: 


G i 


9n 

9 12 

9 13 
<714 
9\ 5 


9n = 0 

912 — ^11 = + iy T 12 + £z T 13) 

<7l3 — ^12 = «/ -1 (£* r 21 + ty T 22 + iz T 2Z) 

914 = ^13 = «7 _1 (£* T 31 + 32 + iz^Zz) 

015 = UT U + VT 12 + WT 1Z - */ _1 (£b91 + ( y?2 + 6?3) 
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(76) 



( 77 ) 



<7 21 
922 

923 , 

9 24 
9 25 


921 — 0 

922 = hi = J^iVxTll + rj v T 12 + 7 h T ls) 

923 = h2 = J^iVxT 21 + f}y^~22 + 7?zT 23 ) 

92 4 = ^23 = J _1 (»?xT31 + 7 ? v r 3 2 + 7 ^X 33 ) 

#25 = UT 2 1 + Vf 22 + Wf 23 - J^iVxqi + »7v?2 + *Ms) 



<731 

932 

933 , 

9 34 
<735 


<731 = 0 

<732 = T31 = «7 _1 (CxTL 1 + Cy r 12 + Cx T 13) 

933 ~ h2 — *^ _1 (Cx r 21 + Cy r 22 + C* t 23) (78) 

<734 = ^33 = «7 _ 1 (C*^31 + Cy r 32 + Cz T 33 ) 

535 = UTzi + VT 3 2 + wf 33 - «7 _1 (C*5l + Cy?2 + C*^) 


9l = 


ar 

a®’ 


ar 


ar 


92 = -a-, * = -a^, 


a = 


l (JL + J*) 

(7 — 1)-Re ' Pr Prt 


From Eqs. (24) and (25), the near-wall shear stress can be written as 

nJk-tLl 

T — T x oa H — PwallUr TT , — r> 

U+ Re y 


(79) 


(80) 


where fi t is an effective turbulent viscosity connecting the wall and the first grid point 


_ y+MiUc 

^ ~ UU+ 

An advantage of Eq. (80) in calculating separated flows is worthy of note: the direction of 
the wall shear stress T wa u is determined by that of the flow velocity U while T wa u calculated 
from Eq. (24) or (2) can only have a positive sign. 

From Eq. (80), the general form of shear stress at the first grid point can be expressed 
as 



t = hail = ATT t (82) 

where 

( fi t !{R e An) if y + > 11.6 

A = < (83) 

( f*waii/{Rt An) otherwise 

U t and An are the tangential component of the resultant velocity and the normal distance 
from the wall, respectively. Eqs. (82) and (83) simply treat the near-wall region as a 
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laminar sublayer (y + < 11.6) and a fully turbulent layer (y + > 11.6). This treatment 
prevents the wall function procedure from producing abnormal results when y + tends to 
zero, such as in the vicinity of separation or reattachment points. 

Similarly, the total heat flux at the first grid point can be written as 

q=-a(T-T wall ) (84) 

where 


a = 


A/[( 7 -l)P, t ] if 2/+ > 11.6 
A/ [(7 — l)P r ] otherwise 
and the heat flux at the wall can be calculated by 


(85) 


qwaii = —a(T - T wa ii ) - Ut ■ T wa u ( 86 ) 

Consider the wall of i7=constant. In this case, only the viscous flux vector G 2 in 
Eq. (75) needs to be modified with the use of the wall functions. In the NPARC code, it 
is calculated by 


dg- 


22 


d 77 

d<7 23 

d 7 ) 
dgiA 
drj 
dff 2S 
dq 


= 922, n — 922,3 

— 922, n ~ 922,a 

= <724, n — 924,b 
= <725, n ~ <725,* 


(87) 


where the subscripts s and n refer to the south and north cell faces at j- 1/2 and j+1/2. 
From Eq. (77), the components of the vector G 2 can be written as 


_ — T— 1/— 2 I _2 , _2\ / X tot^n 

922 J (Vx+Vy + Vz) R ' dl] 


+ 


922 = J 1 (Vx +% + Vz) 


2% 9-tot 


dv 


R e dl] 

2 sf l tot dw 


( 88 ) 




dT 


g 25 = ut 21 + vt 22 + wt 23 + J 1 (t ]1 + % + Vl) a -Q ~ + 


with 


fitot = fl + fit 
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From geometrical consideration, we have 


J \vl •+ -vl + vl) = 


J- 2 {% 


x + vl + Vz) 


J - 1 


AS 2 _ AS 
A ft An 


(Aft = ASAn) (89) 


where Aft and AS are the volume and face area of the control volume, respectively. 

In the wall function approach, all the stresses acting on the cell face considered are 
replaced by the wall shear stress given by Eq. (82). Therefore, for the south wall, Eqs. (88) 
are replaced by 


922 ,. = AS-^u,., = ASAn,., 
922 ,. = AS-^»,. V = ASAu,., 
924. ■ = AS ^ u *,. = ASAn,., 


(90) 


925,1 = ug 2 2,i + v<723,» + wg24,» + ASa(T — T wa a) 

where u t , x , u t , y , u t>z are the x-, y-, z-component of the tangential velocity U t , respectively. 
If n x , n v and n z are the Cartesian components of the unit normal vector at the wall, u t>x , 
u t, y , «t,z can be calculated by 


u t , x = u- n x V n 

Ut,y = V tlyVf i (91) 

u t , z = w- n z V n 

with 


V n = n x u + riyV + n z w 
Similarly, for the north wall, we have 

922, n = —AS ~ — ASAu t>x 

923, n — — AS^^-U t>y = —ASXu t ,y 

< 724 , „ = = -ASAn M 

925, n = Ug 2 2,n + Vg 2 Z,n + W024,n ~ ASa(T - T wall ) 

Walls in the other two directions can be treated in the same way. 


( 92 ) 


18 



The sequence in which the above equations are solved together with the Navier-Stokes 
and turbulence equations in the code is as follows: 

a. Initialize all field values. 

b. Calculate r wall and y + using Eqs. (82) and (25). 

c. Fix the values of K and e at the first grid points using Eqs. (29) and (30). 

d. Solve the turbulence equations. 

e. Calculate q wa u using Eq. (86). 

f. Calculate U e and U+ using Eqs. (26) and (24). 

g. Update fi t using Eq. (81). 

h. Update a using Eq. (85). 

i. Update 522 , 023,924 and 025 using Eq. (90) or (92). 

j. Solve the Navier-Stokes equations. 

k. Return to step b with updated field values. 

The sequence of steps b to k are repeated until the calculation converges. 

Note that by definition, the turbulent eddy- viscosity fi t is zero at the wall, such as in 
the case of the low Reynolds number turbulence models. When using the wall functions, 
Eq. (81) introduces the effective turbulent viscosity which is defined at the wall in the 
turbulence module. Therefore in post-processing, the wall friction coefficient Cf can be 
calculated in the same way as for laminar flows, except that the molecular viscosity fi is 
replaced by the turbulent viscosity fi t at the wall. 
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4 Module Usage 

The present turbulence module (version 2.0) is written based on the NPAEC version 2.2. 
The following are those parts of the code which may require user’s attention when using 
it. 


4.1 Module 

To facilitate identification, all the subroutine names in the module start with CM. In order 
to use the module, the user only needs to call its three subroutines: CM AO, CM ALL and 
CMVRHS, in the NPARC code. 

• Subroutine CMAO. This subroutine transfers from NPARC to the module the param- 
eters to define flow, geometric and boundary conditions. In addition, it has the following 
user-specified parameters: 

JWFAV, KWFAV, LWFAV — Number of grid points above J-, K-, and L-walls at which 
the artificial viscosity is turned off when using the wall function approach. These numbers 
should cover the near- wall region within y + ~ 1000. Currently, they are set to 

JWFAV=5, KWFAV=5, LWFAV=5 

FDEFER — Blending factor in the convection scheme. Its value may vary from 0 to 1 
with the limiting value 0 for the first-order accurate upwind and 1 for the second-order 
accurate HLPA scheme. The solution tends to be more stable, but also more diffusive 
when this factor is reduced. 

BDMAX(i), BDMIN(i) — Upper and lower bounds for the values of K (i=l), e (i=2) and 
fit (i=3). These bounds are introduced for numerical purposes only, that is, to prevent the 
corresponding turbulence quantities from becoming negative or abnormally large during 
the solution process. Currently, they are set to 

BDMAX(l)=1.0E+6, BDMIN(l)=1.0E-8 
BDMAX(2)=1.0E+6, BDMIN(2)=1.0E-8 
BDMAX(3)=1.0E+4, BDMIN(3)=1.0E-3 

which should cover a wide range of the physically meaningful values of K , e and fi t . It is 
to be noted that these values are only valid for the non-dimensional turbulence quantities, 
as defined in the NPARC code. 

• Subroutine CMALL. This is the main subroutine to control the solution sequence 
in the module. The array variable VIST is the turbulent viscosity fi t which is needed in 
NPARC for calculating turbulent flows. The array variables TE, ED and YPS are K, e 
and y + , respectively, which can be used for post-processing. Normally, there is no need 
for user to change this subroutine. 

• Subroutine CMVRHS. This subroutine which is the counterpart of the subroutine 
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VISRHS in NPARC is for introducing the wail function modifications into the right-hand 
side viscous flux terms. There is no need for user to change this subroutine. 


4.2 NPARC 

The authors have made all the modifications necessary for NPARC to use the module. 
The following shows where these modifications are in NPARC. All the alterations are 
marked between C<< and C>> in the code. 

• Namelist TURBIN. The integer parameter IMUTR2 is used to select the turbulence 
models in the module with 

IMUTR2=101 Chien model 

102 Shih-Lumley model 

103 CMOTT model 

A new integer parameter MWALF is introduced to select the near-wall approach with 

MWALF=0 low Reynolds number approach 
1 wall function approach 

Correspondingly, a new statement is added 
in the include file NPARC .INC: 

COMMON /CMOTT /MWALF 
and in the subroutine TURBIN: 

CALL NLGETI(’MWALF\ MWALF). 

• Subroutines FILT1, FILT2, FILTER. An array FAV01 has been introduced into 
each of these subroutines to eliminate the artificial viscosity in the near- wall region when 
using the wall function approach. 

• Subroutine DEPOIN. For using the module, an additional array whose name is 
IPWCM1 is added to IPWLST, and MEMSIZ in NPARC.INC must be increased corre- 
spondingly. 

• Subroutines INITIA, WREST. The model identifier ITURB for each K-e turbu- 
lence model in the module is given an integer value greater than 100. To reflect this 
expanded choice for turbulence models, the read and write statements for the turbulent 
quantities in these two subroutines are modified as 

IF(ITURB.EQ.4 .OR. 

& ITURB.EQ.5 .OR. 

& ITURB .EQ. 7 .OR. 

& ITURB. GT. 100) 

& READ(2) or WRITE(4) (values of K, e, n t ) 
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• Subroutine MUTURB. Here, the subroutines CMAO and CM ALL of the module 
are called as follows: 


ELSE IF(ITURB.EQ.7) THEN 

CALL TKWSTEP(BIGA(IPWR),BIGA(IPWAK),BIGA(IPWEPS) 
&,BIGA(IPWS1)) 

C« 

ELSE IF(ITURB.EQ.101 .OR. ITURB.EQ.102 .OR. ITURB.EQ.103) THEN 
CALL CMAO(JMAX,KMAX,LMAX,NMAX,NTURB,NRLX,NSPRT,ITURB 
&,MWALF 

&, RE, C2B, GAMMA, RPR, RPRT,DT,DTCAP,IVARDT 
& ,TUIN 1 ,TUIN 2 ,TUIN 3 ,TMUIN 1 ,TMUIN2 ,TMUIN3 
&,NJPAT,JPJ2,JPJM,JPK2,JPKM,JPL2,JPLM 
&,NKPAT,KPJ2,KPJM,KPK2,KPKM,KPL2,KPLM 
&,NLPAT,LPJ2,LPJM,LPK2,LPKM,LPL2,LPLM 

&,NJSEG,JLINE,JKLOW,JKHIGH,JLLOW,JLHIGH,JTYPE,JSIGN,JEDGE 

&,TEMPJ 

&,NKSEG, KLINE, KJLOW,KJHIGH,KLLOW,KLHIGH,KTYPE,KSIGN, HEDGE 
&,TEMPK 

&,NLSEG,LLINE,LJLOW,LJHIGH,LKLOW,LKHIGH,LTYPE,LSIGN, LEDGE 
&,TEMPL) 

C 

IPWR3=IPWS1+2*NXYZ 

IPWR4=IPWR3+NXYZ 

CALL CMALL(BIGA(IPWX ),BIGA(IPWY ),BIGA(IPWZ ) 

&,BIGA(IPWXX ),BIGA(IPWXY ),BIGA(IPWXZ ) 

&,BIGA(IPWYX ),BIGA(IPWYY ),BIGA(IPWYZ ) 

&,BIGA(IPWZX ),BIGA(IPWZY ),BIGA(IPWZZ ),BIGA(IPWQ ) 
&,BIGA(IPWR ),BIGA(IPWRU ),BIGA(IPWRV ),BIGA(IPWRW ) 
&,BIGA(IPWE ),BIGA(IPWAK ),BIGA(IPWEPS),BIGA(IPWTMU) 
&,BIGA(IPWR3 ),BIGA(IPW28 ),BIGA(IPW29 ),BIGA(IPWR4 ) 
&,BIGA(IPWVDT),BIGA(IPWS1 ) 

&,BIGA(IPWS1 + 4*NXYZ ),BIGA(IPWSPT),BIGA(IPWCOF) 
&,NC,SUMN,TMAX) 

C» 

ELSE 


END IF 


• Subroutine STPF3D. The subroutine CMVRHS of the module is called here. 
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5 Demonstration Examples 

Two examples are given here to demenstrate how to use the NPARC code with the 
turbulence module. Other application examples can be found in Yang et al. (1995). 


5.1 Flat Plate 

Turbulent boundary layer flow over a flat plate with zero pressure gradient was selected 
as the first test case for code validation. The solutions of both the low Reynolds number 
(LR) and the wall function (WF) approaches were compared with the experimental data 
(Exp) of Wieghardt and Tillmann (1951). Fig.l shows the flow geometry and boundary 
conditions used in the calculation. For the wall function approach, gnd points were 
111x55x5, the first grid points above the wall had the y + value of 60 and 14 grid points 
in the x-direction were located before the leading edge of the flat plate. For the low 
Reynolds number approach, grid points were 111x81x5 with the distribution of x-points 
being the same as in the wall function case and the first y+ being 0.3. Since the NPARC 
code is for compressible flows while the experiment to be compared was for incompressible 
flows, a freest re am Mach number of 0.2 was chosen. All calculations started from an initial 
field with U = 0.2, V = 0,/x t = l,K = 0.005 and e = 0.09 R*K 2 /fi t . The detailed NPARC 
input data are given in Appendix 2. Fig.2 shows the convergence history of the Chien 
model. The other two models had the similar convergence behavior. For comparison, 
the result of the Raldwin-Lomax model in the NPARC code is also given in Fig-2. The 
low Reynolds n um ber Chien model had a rather poor convergence behavior: its L2 norm 
of residual was only reduced to 4.7 xl0~ 8 after 20000 iterations, and further down to 
8.9 xlO~ 10 after another 60000 iterations. However, the solutions at both the 20000th 
and 80000th iterations were almost identical. Fig.3 shows the wall friction coefficient C f 
versus the Reynolds number based on the momentum thickness of boundary layer Re$. 
It can be seen that for the the low Reynolds number approach, the results of both the 
Shih-Lumley and CMOTT models are almost the same and slightly better than that of 
the Chien model, while for the wall function approach, the results of both the Chien and 
Shih-Lumley models are almost identical and slightly better than that of the CMOTT 
model. Regarding the computational cost of the Chien model, 1000 iterations on the Cray 
YMP computer took 1566 seconds for the wall function approach and 2320 seconds for 
the low Reynolds number approach. 

5.2 Transonic Diffuser 

This case was taken from the experiment of Salmon et al. (1983). Fig.4 shows the 
computational domain of the diffuser. In the experiment, different measurements were 
taken, ranging from no-shock to strong-shock conditions. During the code validation, we 
calculated three cases with no-shock, weak-shock and strong-shock, respectively. Here, 
we only show the strong-shock case which is most difficult to calculate accurately. All 
calculations started from an initial field with U = 0.44, V = 0 ,p = 0.5143, fit = 10 , if — 
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0.0001 and e = 0.09R e K 2 / p t . For the low Reynolds number approach, the grid had 
81x81x5 points and its x-point distribution was the same as that used by Dudek (1995), 
and for the wall function approach, the grid had 81x51x5, differing Rom the former one 
only in the y direction. The detailed NPARC input data are given in Appendix 2. Fig.5 
shows the convergence history, and it can be seen that the wail function approach has 
a better convergence behavior than the low Reynolds number approach. The pressure 
distributions are shown in Fig. 6 for the top wall and in Fig.7 for the bottom wall. For 
the low Reynolds number approach, the CMOTT model gave the best results, especially 
for the location of the shock wave and the pressure recovery after the shock wave, while 
for the wall function approach, the results of the CMOTT model were slightly worse than 
those of the other two K-e model. Fig.8 shows the streamwise mean velocity profile at the 
section x/Xp e /=4.6, and it can be seen that the results of all the three K-e models with 
the wall functions are nearly the same, and better than those of the low Reynolds number 
counterpart. For the computational cost of the Chien model, 1000 iterations took 1088 
seconds for the wall function approach and 1720 seconds for the low Reynolds number 
approach. 
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Fig.l Flow geometry and boundary conditions 



Fig.2 Convergence history 
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Fig.3 Wall friction coefficient 
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Fig.4 Computational domain of diifuser 



Fig. 5 Convergence history 
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Fig.7 Pressure distribution at the bottom wall 
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7 APPENDIX 2: NPARC Input Files 

• Flat Plate Case, Wall Function Approach 

$ INPUTS 


pref 

= 

2116.80 

nblock 

= 

1 

tref r 

= 

519.0 

nmax 

= 

1000 

vrat 

= 

1 . 000001e+00 , 

nc 

= 

-1 

tsuth 

= 

0 . 198680e+03 , 

nsprt 

= 

50 

re 

= 

5 . 000000e+06 , 

np 

= 

960000 

pr 

= 

0 . 720000e+00 , 

ifxprt 

= 

0 




if xplt 

= 

0 

dtcap 

= 

5.0 

12plot 

= 

-1 

pcqmax 

= 

1.000000e+01 , 

iplot 

= 

0 

splend 

= 

0.500000e+00 , 




stopl2 

= 

1 . 000000e-15 , 

numdt 

= 

0 

stoptr 

= 

100.0 

ivardt 

= 

2 

xmach 

= 

0 . 200000e+00 , 

ifiltr 

= 

2 




ispect 

= 

1 

gamma 

= 

1 . 400000e+00 , 

lrest 

= 

1 

dis2 

= 

0.250000e+00 , 

mbord 

= 

1 

dis4 

= 

0 . 640000E+00 , 

ibord 


1 


$END 

$TURBIN 

imutur = 101, 

imutr2 = 101 , 

mwalf = 1, 

nturb = 0 , 

order = 1 . 0 , 

nrlx = 1, 

ifmax = 3, 

prt = 0.9, 

tuinl = 0 . 02 , 
tmuinl = 10 . 0 , 

$END 

5BL0CK 



invisc(l) = 

1, 

invisc(2) = 

i, 


invisc(3) = 1, 




lamin(l) 

0, 

lamin(2) = 



lamin(3) = 0, 




npseg = 

1 , 

nbcseg 

= 

7, 





$END 










$PRTSEG 










jklpi (1, 1, 1) 

= 

1, 55, 

55, 







jklpi (1,2,1) 

= 

1 , 55, 

10, 







jklpi (1,3,1) 

= 

3, 3, 

1, 







ipord (1,1) 

= 

2, 1, 

3, 






$END 










1 

1 2 

54 

1 5 

0 


1 

0.7345000 

1.0080000 

-1 

111 

111 2 

54 

1 5 

0 

- 

■1 

0.7143000 

1.0000000 

0 

1 

14 1 

1 

1 5 

50 


1 




15 

111 1 

1 

1 5 

60 


1 



50 

1 

111 55 

55 

1 5 

50 

- 

•1 




1 

111 1 

55 

1 1 

50 


1 




1 

111 1 

55 

5 5 

50 

- 

’1 
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• Flat Plate Case, Low Reynolds Number Approach 

$ INPUTS 


pref 

= 

2116.80 

r 

nblock 

i , 

trefr 

= 

519.0 

/ 

nmax 

= 1000 , 

vrat 

= 

1 . 000001e+00 

/ 

nc 

-1 , 

tsuth 

= 

0 . 198 680e+03 

/ 

nsprt 

50 , 

re 


5 . 000000e+06 

r 

np 

= 960000, 

P r 

= 

0 . 720000e+00 

t 

ifxprt 

0 , 





ifxplt 

0 , 

dtcap 

= 

1.0 

r 

12plot 

-1 , 

pcqmax 

= 

1 . 000000e+01 

f 

iplot 

0 , 

splend 

= 

0 . 500000e+00 

/ 



stopl2 

= 

1 . 000000e-15 

/ 

numdt 

0 , 

stoptr 

= 

100.0 

/ 

ivardt 

2 , 

xmach 

= 

0 . 200000e+00 

f 

if iltr 

2 , 





ispect 

1 , 

gamma 

= 

1 . 400000e+00 

r 

lrest 

1 , 

dis2 

= 

0 . 250000e+00 

t 

mbord 

1 , 

dis4 

= 

0 . 640000E+00 

t 

ibord 

1 , 


$END 

$TURBIN 

imutur = 2, 

imutr2 - 101 , 

mwalf = 0, 

nturb = 4001, 
order = 1.0, 

nrlx = 100, 

ifmax = 3, 

prt = 0.9, 

tuinl = 0.02, 
tmuinl = 10.0, 

$END 

$BLOCK 

invisc(l) = 1, invisc(2) = 1, invisc(3) = 1, 

lamin(l) = 0, lamin(2) = 1, lamin(3) = 0, 

npseg = 1, nbcseg = 7, 

$END 

$PRTSEG 

jklpi(l,l,l) = 1, 81, 81, 

jklpi (1,2, 1 ) = 1, 81, 10, 

jklpi (1,3,1) - 3, 3, 1, 

ipord(l,l) = 2, 1, 3, 

$END 

1 1 2 80 1 5 0 1 0.7345000 1.0080000 -1 

111 111 2 80 1 50 -1 0.7143000 1.0000000 0 

1 14 1 1 1 5 50 1 

15 111 1 1 1 5 60 1 70 

1 111 81 81 1 5 50 -1 

1 111 1 81 1 1 50 1 

1 111 1 81 5 5 50 -1 
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• Transonic Diffuser Case, Wall Function Approach 


$ INPUTS 


pref 

= 

2819.5 

r 

nblock 

= 

1 


trefr 


525.6 

/ 

nmax 

= 

1000 


vrat 

= 

0.100000e+01 

/ 

nc 

= 

-1 


tsuth 


0 . 198680e+03 

/ 

nsprt 

= 

50 


re 

= 

1.339171e+06 

/ 

np 

= 

900000 


pr 

= 

0.720000e+00 

/ 

ifxprt 

= 

0 






ifxplt 

= 

0 


dtcap 

— 

0.5 

/ 

12plot 

= 

-1 


pcqmax 

= 

10.0 

r 

iplot 

= 

0 


splend 

= 

0.50 

r 





stopl2 

= 

1 . 000000e-15 

/ 

numdt 

= 

0 


stoptr 

= 

100.0 

t 

ivardt 

= 

2 


xmach 

= 

0.10 

t 

if iltr 

= 

2 






ispect 

= 

1 


gamma 

= 

1.40 

t 

lrest 

= 

1 


dis2 

= 

0.25 

t 

ibord 

= 

1 


dis4 


0.64 

r 

mbord 

= 

1 


$END 








$TLJRBIN 








imutur 

= 

2, 






imutr2 

= 

101, 






mwalf 


1, 






nturb 

= 

0, 






order 

= 

1.0, 






nrlx 


1/ 






ifmax 


2, 






prt 

= 

0.9, 






tuinl 

= 

0.02, 






tmuinl 

= 

10.0, 







$END 


$ BLOCK 

invisc(l)=l, invisc(2)=l, invisc(3)=l, 
lamin(l) =1, lamin(2) =1, lamin(3) =1, 
npseg =1, nbcseg =6, 

SEND 

SPRTSEG 



jklpi ( 1 , 1, 

1) = 

1, 

41 , 

41, 





jklpi (1,2, 

1) = 

1, 

51, 

10, 





jklpi (1,3, 

1) = 

3, 

3, 

1, 





ipordd, 1) 

= 

2, 

1, 

3, 




$END 









1 

1 2 

50 

1 

5 

0 

1 

0.7143000 1.0000000 

-1 

81 

81 2 

50 

1 

5 

0 

-1 

0.5143000 0.5000000 

0 

1 

81 1 

1 

1 

5 

60 

1 


25 

1 

81 51 

51 

1 

5 

60 

-1 


25 

1 

81 1 

51 

1 

1 

50 

1 



1 

81 1 

51 

5 

5 

50 

-1 
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• Transonic Diffuser Case, Low Reynolds Number Approach 


$ INPUTS 


pref 

= 

2819.5 

/ 

nblock 

= 

1 

r 

trefr 

= 

525.6 

, 

nmax 

= 

1000 

r 

vrat 


0 . 100000e+01 

t 

nc 

= 

-1 

f 

tsuth 

= 

0 . 198680e+03 

f 

nsprt 

= 

50 

r 

re 

= 

1 . 339171e+06 

t 

np 

= 

900000 

/ 

pr 

= 

0 . 720000e+00 

r 

ifxprt 

= 

0 

r 





ifxplt 

= 

0 

r 

dtcap 

= 

0.5 

, 

12plot 

= 

“1 

t 

pcqmax 

= 

10.0 

r 

iplot 

= 

0 

f 

splend 

= 

0.50 

/ 





stopl2 

= 

1 . 000000e-15 

/ 

numdt 

= 

0 

r 

stoptr 

= 

100.0 

/ 

ivardt 

= 

2 

t 

xmach 

= 

0.10 

r 

ifiltr 

= 

2 

f 





ispect 


1 

t 

gamma 

= 

1.40 

t 

lrest 

= 

1 

/ 

dis2 


0.25 

t 

ibord 

= 

1 

r 

dis4 

= 

0.64 

f 

mbord 

= 

1 

/ 

$END 








$TURBIN 








imutur 

= 

2, 






imutr2 


101, 






mwalf 

= 

o, 






nturb 

= 

0, 






order 

= 

1.0, 






nrlx 

= 

1, 






ifmax 

= 

2, 






prt 

= 

0.9, 






tuinl 

= 

o 

o 

ro 







tmuinl = 10.0, 

SEND 

SBLOCK 

invisc(l)=l, invisc(2)=l, invisc(3)=l, 
lamin(l) =1, lamin(2) =1, lamin(3) =1, 
npseg =1, nbcseg =6, 

SEND 

SPRTSEG 



jklpi (1, 1, 

i) 

= 

1/ 

41, 

41, 






jklpi (1,2, 

i) 

= 

1, 

81, 

10, 






jklpi (1, 3, 

i) 


3/ 

3, 

1, 






ipord (1, 1) 


— 

2, 

1, 

3, 





$END 











1 

1 2 


80 

1 

5 

0 

1 

0.7143000 

1.0000000 

-1 

81 

81 2 


80 

1 

5 

0 

-1 

0.5143000 

0.5000000 

0 

1 

81 1 


1 

1 

5 

60 

1 



40 

1 

81 81 


81 

1 

5 

60 

-1 



41 

1 

81 1 


81 

1 

1 

50 

1 




1 

81 1 


81 

5 

5 

50 

-1 
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